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Abstract 

Event-Driven Particle Dynamics is a fast and precise method to simulate partic- 
" ^ . ulate systems of all scales. These advantages arise from the analytical solution 

of the dynamics required by the discrete-potential models used. Despite the 
high precision solution, the finite calculation-precision of computers will still 
cause the simulation to enter invalid states which, if left unchecked, can lead 
to unresolvable errors. In this work, the treatment of these marginal invalid- 
states is discussed and a general event-detection algorithm is proposed which 
stably handles these situations. This requires a definition of the dynamics of in- 
valid states and leads to improved algorithms for event-detection in spherically 
symmetric systems, including the well-established hard-sphere and square-well 
models. Finally, the Event-Driven Particle Dynamics technique is extended to 
allow the study of systems with complex spherical-mesh boundary conditions 
and distance constraints as a demonstration of the generality of the proposed 
algorithm. 

Keywords: DEM; event driven; molecular dynamics; square well; hard sphere; 
stepped potentials; 

Event-driven particle-dynamics (EDPD) algorithms are the oldest particle 
dynamics technique [1] and have found application in a wide range of fields, from 
predicting vapour- liquid equilibria [2] to protein folding studies [3] and the design 
of granular vibration dampers [4]. These applications of EDPD differ only in the 
selection of models used to represent the inter-particle forces. For conservative 
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Figure 1: The potential energy u(r) of the square-well model as a function of separation 
distance, r, between two particles. The inter-particle energy only changes in discrete steps at 
the distances r = a and r = Act. The numbered arrows indicate the four event types of the 
model generated by these steps: (1) Capture (enter the well), (2) Release (escape the well) 
and Bounce (remain in the well), and (3) Core [8]. 

inter-particle forces, these are often specified through an inter-particle energetic 
potential. The EDPD technique can only be applied to models where this poten- 
tial changes at discrete separation distances; for example, see the inter-particle 
potential of square- well molecules in Fig. 1. Despite this restriction, the EDPD 
technique remains a general approach to particle simulations as discontinuous 
potentials can be "stepped" to accurately approximate more conventional model 
systems, such as Lennard-Joncsium[5, 6], or directly fit to physical data[7]. 

A significant advantage of discontinuous particle models is that they are 
amenable to theoretical analysis [8] and are at the heart of the most successful 
molecular theories, including thermodynamic perturbation theory[9] and kinetic 
theory [10]. The second major advantage of using discrete potential models is 
that the EDPD simulation technique is fast and does not suffer from truncation 
error due to the analytical nature of the dynamics [11]. Despite these significant 
advantages, the EDPD method is not as popular as the ubiquitous time-driven 
techniques [12]. This is due, in part, to the intricate and complex nature of 
the EDPD algorithm and its high sensitivity to numerical and implementation 
errors. To illustrate these points, the general time-stepping algorithm is briefly 
outlined and contrasted against the event-driven approach. 

In time-stepping particle dynamics algorithms, the inter-particle potential, 
u(r), is used to define the inter-particle forces. The sum of all of the forces 
acting on a particle is then inserted into Newton's equation of motion and nu- 
merically integrated to solve for the trajectory of the particle. This approach 
is straightforward but it requires a careful selection the numerical integration 
technique and corresponding time step [11]. Time-stepping simulations must 
be performed with continuous potentials as it becomes impossible to define a 
force suitable for numerical integration when a discontinuous model, such as the 
square- well model in Fig. 1, is used. This is a result of the interparticle force 
being equal to the negative gradient of the potential, which is not easily defined 
at a discontinuity. Instead, whenever a pair of particles transition a disconti- 
nuity in their potential, the interaction dynamics must be treated separately 
to ensure momentum and energy are conserved. Fortunately, as the potential 



2 



is flat between discontinuities, there are no inter-particle forces acting between 
transitions and Newton's equation of motion can be solved analytically for all 
particles. This analytical dynamics allows the calculation of the times that par- 
ticle pairs will transitions across discontinuities, known as "events" , without 
calculating intermediate states or using numerical integration. The simulation 
can then use the analytical dynamics to skip forward to the time of the next 
event and separately process the interaction at the exact time the event occurs. 
The simulation algorithm has become "event-driven" and the calculation time 
of the simulation is now primarily spent on detecting and executing events [13] 
rather than numerical integration. 

Although the analytical nature of event-driven dynamics results in a high 
precision solution of the equations of motion, a simulation will slightly under or 
over estimate the time of the transition due to round-off errors in the calcula- 
tions. This may cause the simulation to move particles beyond discontinuities 
which are energetically inaccessible, leading to the particle pair entering an 
invalid state. This problem is exacerbated in particle models including rota- 
tional degrees of freedom as numerical root finding techniques must be used to 
detect events[14, 15, 16] and these may miss glancing interactions, again lead- 
ing the system to access invalid states. There are some special cases for which 
guaranteed numerical event-detection routines are available, such as those avail- 
able for infinitely-thin hard rods[14], but these techniques are not available in 
general[15] and are still vulnerable to round-off errors. Although these round-off 
errors are small in magnitude, once particles enter invalid states the dynamics 
of the model is no longer defined. As simulations will accidentally access invalid 
states, it is imperative that event-detection algorithms carefully define the dy- 
namics in these regions to ensure that these invalid states will resolve in such 
a way as to minimise their effects. It is interesting to note that this ambigu- 
ity in the overlapped dynamics has also led to some difficulties in theoretical 
treatments in the past [17], but these were resolved by carefully ensuring that 
theoretically inaccessible states have no bearing on the final results. 

Despite EDPD's long history, it is difficult to find event-detection algorithms 
which carefully treat invalid states for even the simplest particle models. Pub- 
lished algorithms will often implicitly rewind simulation time until the invalid 
state is eliminated[l, 11]. This approach is not stable as the rewind is unchecked 
and may result in new invalid states appearing for other particle pairs. Careful 
implementations have suggested disabling interactions for particles in invalid 
states [18], but this fails to prevent even minor precision errors from deterio- 
rating into unrcsolvable errors. An excellent discussion on these difficulties is 
presented by Poschel and Schwager[19], who propose to combat round-off errors 
by biasing the movement of particles so that detected interactions occur at a 
small distance from the invalid state, protecting against round-off errors. This 
approach is unacceptable as it modifies the dynamics of the model in the valid 
states, it will not prevent invalid states developing in confined geometries, and it 
is unable to resolve existing minor invalid states. It is clear that although these 
minor errors are well-known and are of a small magnitude the current published 
methods do not handle these errors in a satisfactory way. 
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It is crucial that EDPD event-detection routines are carefully constructed 
to ensure the algorithms are stable with regard to invalid states and always 
minimise their duration and effect on the dynamics. In this work, a general 
methodology for creating stable event-detection routines for use in event-driven 
simulations is developed. In the following section, the general approach is de- 
rived and its justification is presented through a discussion of spherically sym- 
metric potentials. In Sec. 2, the solutions for spherically symmetric potentials 
are derived. This results in improved expressions for common EDPD systems, 
as well as a stable set of algorithms for systems including gravity. Finally, in 
Sec. 3 the results of the paper are summarised and a discussion of future work 
is presented. 

1. General Methodology 

In this section, the general approach for generating stable event-detection 
algorithms is outlined. This discussion is motivated through a consideration of 
the square- well model, as illustrated in Fig. 1. This model is ideally suited for 
this discussion as it is relatively simple yet it contains all of the event-detection 
difficulties of all spherically-symmetric discrete potentials. Although the ideas 
outlined here readily generalise to asymmetric systems, a demonstration of this 
will be reserved for a future publication. 

When searching for the next event to occur in a square- well system, all pairs 
of particles can be considered separately for events as there is a vanishing proba- 
bility of interactions between three or more particles [1 1] . All pairings of particles 
are tested for events and the earliest event is, by definition, the next event to 
occur in the system. This simplifies the event-detection for the system into the 
detection of events for pairs of particles. Considering a single pair of particles, 
every discontinuity in the inter-particle potential can be individually tested for 
transitions and, again, the earliest transition is the next candidate event for 
the pair. Thus the problem of event-detection for the whole system of particles 
reduces to a search for an algorithm for detecting transitions across a single 
discontinuity for a single pair of particles. For spherically symmetric potentials, 
this problem simplifies further to finding the earliest time a pair of particles will 
reach a separation distance equal to the location of the discontinuity. 

It is crucial to note that only the discontinuities which are accessible to the 
pair of particles need to be tested for transitions. Which discontinuities are 
accessible depends on the state of the particle pair. For example, uncapturcd 
square-well particles are outside of the well of the potential and only need to 
be tested for Capture events, marked 1 in Fig. 1. Once a particle pair becomes 
"captured", the pair has entered the well of the potential and only needs to 
be tested for Core and Bounce/Release events, marked 2 and 3 in Fig. I re- 
spectively. The captured/uncaptured state of the particle pair also alters which 
states of the model are defined as invalid, as illustrated in Fig. 2. Uncapturcd 
particles nominally lie outside the bounds of the well (see Fig. 2a) but captured 
particle pairs nominally lie inside the well (see Fig. 2b and c). Considering one 
discontinuity at a time, there are two classes of event-detection algorithm which 
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Figure 2: An illustration of the invalid states and interactions for a pair of (a) uncapturcd 
and (b— c) captured square- well particles. The coordinate frame is centred on particle j, and 
the trajectory of particle i is indicated with a dotted line. The next interaction of the particle 
pair is marked using a filled circle. The interaction occurs at the distance of the discontinuity, 
denoted cy, which is listed underneath each plot. The invalid states of the particle pair are 
shaded grey. 



may be distinguished: The first is where the valid state lies on the positive side 
of the discontinuity and the second is where the valid state lies on the nega- 
tive side. For example, the Capture event is valid on the positive side of the 
discontinuity at Acr + (Fig. 2a) but the Bounce/Release event is valid on the 
negative side of the discontinuity at Aer~ (Fig. 2b). Now that the problem of 
event-detection is well-defined, and the invalid states of the discontinuity can be 
identified from the sign of the discontinuity, an event-detection algorithm which 
is stable with respect to invalid states can be derived. 

To broaden the application of the results obtained here, each particle of the 
system is considered to be under a constant acceleration. Particle systems which 
include constant accelerations, such as gravity, are some of the most numerically 
sensitive due to the frequent re-collisions of particles and the nonlinear particle 
trajectories [18]. It is also common in time-stepping simulations to use particle 
meshes, which are an array of static particles with infinite mass and zero accel- 
eration, to create the boundaries of the simulation[19]. Thus the acceleration 
constant is also allowed to vary in value for each particle to enable these meshes 
to be implemented in EDPD. Considering the motion of a single particle i be- 
tween events, the equation of motion for the particle can be solved analytically 
to yield 

At 2 

r l (t + At) = r l {t) + Atv l {t) + —a i (1) 

where T"i(f) is particle i's position at the current time t, Vi(t) is it's current 
velocity, and is it's constant acceleration. Equation (1) is used in simulations 
to move particles forward in time; however, to detect events between a pair of 
particles, it is convenient to transform the coordinate system relative to the 
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second particle, j. The equation of motion for particle i then becomes 

At 2 

m (t + At) = m (t) + At Vij (t) + —a^ (2) 

where = r,— Tj is the position of particle i relative to particle j, Vij — Vi—Vj 
is the relative velocity, and a,j = — a,j is the relative acceleration. 

Assuming there is a discontinuity in the potential at a separation distance of 
cry , the detection of the next transition of the particles across this discontinuity 
can be expressed as a search for the smallest positive value of At which is a 
solution to the following equation: 

r ij (t + At)=a ij r ij (t + At) (3) 

where f%j is the unit vector in the direction of the separation vector . Finding 
a solution to Eq. (3) in its current form is difficult, but it can be greatly simplified 
by squaring both sides of the equation to yield: 

(r y -(t + At)) 2 -4=0. (4) 

The search for transition times across a discontinuity has now been transformed 
into a search for the roots of the left-hand side of Eq. (4). Isolating these terms, 
we can define the overlap function, /, for spherically symmetric potentials as 
follows 

/+(* + At) = -f-(t + At) = ( Py (t + At)) 2 - 4, (5) 

where the superscript on the overlap function is used to denote the two pos- 
sible choices for the sign of the overlap function. In this work it is required 
that the definition of the overlap function is chosen carefully so that the func- 
tion is negative for particle pairs in an invalid state, and positive in all valid 
states. Therefore, the positive definition, / + , is used to detect interactions on 
the positive side of the discontinuity, such as Capture {a^ — A a + ) and Core 
(<7ij = cr + ) events. The negative definition, /~, is used to detect interactions 
on the negative side of the discontinuity, such as the Bounce/Release events 
(cry = Aer~). The use of an overlap function is common in studies on event- 
driven asymmetric-potential systems, where it has also been called the overlap 
potential[15] or the indicator function[16]. 

The key motivation for defining an overlap function is to transform the 
search for the transition times of the particle pair across a discontinuity at <Xy 
into a search for the roots of the overlap function, /(t). In addition, the sign 
of the overlap function can be used as a test if the particle pair is in an invalid 
state. On inspection of Eq. (5), it is also clear that the overlap function is also 
a measure of the distance of a pair of particles from the discontinuity at cry. 
Crucially, these properties allow the derivative of the overlap function to be 
used as a test if an invalid state (/(t) < 0) is approaching a valid state with 
time (/(t) > 0) or not (/(f) < 0). Thus, a stable EDPD algorithm may be 
defined as one which ensures the overlap function is increasing for all particles 
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in invalid states. This definition of a stable algorithm requires it to generate a 
response to cases where the particles are in an invalid state (f(t) < 0) which are 
not improving (f(t) < 0), but what response is appropriate is not immediately 
clear. 

Whenever a pair of particles are about to enter an invalid state, the dynamics 
of the model must require that an interaction occurs and that these interactions 
must cause the system to move towards a valid configuration. For example, 
the repulsive Core interaction (see Fig. 2b) causes a pair of captured square- 
well particles with contacting cores to begin to move apart, away from the 
invalid state of overlapping cores. The Bounce interaction prevents particles 
escaping the well until they have sufficient energy for a Release event. Therefore, 
interactions are, by their very nature, stabilising actions and here it is proposed 
that they may be used to attempt to stabilise pairs of particles in invalid states 
where the overlap function is decreasing. A stable event-detection algorithm 
can then be obtained by applying the following rule: 

Given an overlap function, f, for an accessible discontinuity, a pair of par- 
ticles will interact if, at any time t, the particles are overlapped or in contact, 
fit) < 0, and the overlap function is decreasing, f(t) < 0. 

This algorithm is fundamentally stable as it prevents invalid states from 
deteriorating at any time by executing an interaction. At the same time, this 
rule will also not allow the overlap functions for any other particle pair in the 
system to become negative, or decrease while negative. A general implementa- 
tion of this rule is given in Algorithm 1. Although this implementation appears 
straightforward, the FindNextRoot function may be extremely complex[15]. 

1 if /(*) > then 

2 | return FindNextRoot (/, t); 

3 else 

if fit) < then return t; 

tderiv < — FindNextRoot(/, i); 
if tderiv = 00 then return oo; 
if f (tderiv) > then 

return FindNextRoot (/, tderiv)', 
else 

return tderiv] 
end 
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5 
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10 
11 

12 end 

Algorithm 1: The general stable EDPD algorithm for event detection. This 
algorithm returns the time of the next transition across the discontinuity at 
a separation of o~ij for particle i and j. The FindNcxtRoot(/, t) function 
returns t + At where At is the smallest positive non-zero value which satisfies 
fit + At) = 0. If there is no such root, it returns positive infinity which 
corresponds to no interaction. 

The general algorithm will prevent invalid states from deteriorating; how- 
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ever, the interaction time returned in Algorithm 1 on Line 10 is problematic if 
it is executed by the simulation. This line schedules an interaction at a turning 
point in the overlap function to prevent the overlap function decreasing while 
the system is in an invalid state. As the overlap function is a measure of parti- 
cle separation, this is a time when the relative velocity along the line of contact 
between the particles is zero. Almost all interaction rules require a non-zero rel- 
ative velocity in the direction of contact to have any effect, and so this case may 
result in an infinite number of ineffectual interactions being scheduled in a finite 
time. This case will occur either in poorly initialised systems, or in dissipative 
systems, where it may appear during an inelastic collapse [20]. The general al- 
gorithm is still susceptible to dynamical-arrest effects, such as inelastic collapse; 
however, these effects are a phenomenon of the model and do not arise due to al- 
gorithmic/implementation difficulties. Modifications to the model to overcome 
these difficulties arc not treated here and are discussed elsewhere[20, 21, 22, 23]. 
As these dynamical arrests are only seen in molecular/elastic systems when the 
initial configuration is invalid, it is possible to simply report it as an unrecover- 
able error. If automatic overlap resolution is desired, this particular case may 
be treated by rotating the relative velocities of the particles or by the injection 
of some additional energy into the system, as previously suggested for inelastic 
collapsc[21]; however, this is a modification of the model of the system and it is 
not recommended for "production" simulations. 

Although Algorithm 1 is a general approach to creating stable EDPD al- 
gorithms, it is not optimal for all cases. The FindNextRoot function is also 
particularly complex to define in systems involving gravity or asymmetry as 
the overlap function may be a high order polynomial or transcendental. In the 
following sections, the general algorithm is applied to spherically symmetric sys- 
tems with gravity to demonstrate that optimal implementations of Algorithm 1 
can be carried out in these important cases. 

2. Solution for Spherically Symmetric Interactions 

Spherical potentials are perhaps the most important shapes to be able to 
detect events for, as most EDPD results are obtained on spherically symmetric 
potentials and spheres are also often used to bound more complex shapes to 
improve performance [24] . It is therefore crucial that event detection for spherical 
potentials is implemented in a stable and optimal fashion. If the equation of 
motion for a pair of particles (Eq. (2)) is substituted into the expression for the 
overlap functions (Eq. (5)), the following expressions arc obtained 

At 4 

f+(t + At) = -f~(t + At) = —a% + At 3 a l} • «y 

+ At 2 (a l3 ■ rij + v%) + 2 At Vij ■ m +r%-al (6) 

f + (t + At) = -f~ (t + At) = At 3 a% + 3 At 2 a rj ■ 

+ 2 At (a, 3 • nj + v 2 j) + 2 Vij ■ m (7) 




Figure 3: The possible classes of trajectories for particles with zero relative acceleration 
(djj = 0). The dotted line represents the trajectory of particle i relative to the centre 
of particle j. The solid line marks sections of the trajectory which meet the condition for 
interaction. Roots of the overlap function (Eq. (8)), corresponding to inter-particle separations 
of cry, are marked with circles and roots of the time derivative of the overlap function (Eq. (9)), 
corresponding to maxima and minima of the overlap function, are marked with triangles. 
Filled triangles or circles indicate that the root marks the start of a section of the trajectory 
which satisfies the conditions for an interaction. 

where all variables are evaluated at a time t. The task of detecting events has 
become a task of finding the roots of Eq. (6). Although the roots of quartic 
equations such as Eq. (6) can be solved using radicals, in general there is no 
strategy which is guaranteed to correctly detect all roots [25]. Fortunately, in the 
majority of applications all particles have zero relative accelerations (a.^ = 0). 
This greatly simplifies the problem and is discussed in the following subsection, 
before moving on to the solution for non-zero relative accelerations. 

2.1. Zero relative acceleration (a^ =0) 

The case of zero relative acceleration occurs if there is no acceleration or 
both particles are under the influence of the same acceleration, e.g. gravity. In 
this case the overlap function reduces to a quadratic and the time derivative 
becomes linear in time 

f+ j=0 (t + At) = -/" , =0 (* + At) = At 2 v% + 2 At Vij ■ r y + r% - a% (8) 

/ay=o(* + At ) = -/a„=o(* + At) = 2 AtV% + 2 Vij ■ TV (9) 

Before attempting to derive the interaction algorithms, it is convenient to draw 
all possible roots and interaction configurations for both overlap functions, as in 
Fig. 3. These diagrams help to quickly observe which parts of a trajectory will 
generate interactions and to identify patterns in the roots of the overlap function 
and its derivative. This aids in the optimisation of the general algorithm for 
these specific cases. 

First, the event detection for interactions on the positive side of the discon- 
tinuity are considered. This is the event-detection rule for Core 1 and Capture 
events in the square well model, and it corresponds to the detection of the the 
roots of the positive overlap function, f^ ij=0 (sec Eq. (8) and Fig. 3a-b). There 
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is only one configuration of roots where the particles cross the discontinuity and 
the overlap function has two roots, which is illustrated in Fig. 3a. The gradient 
of the overlap function, / + , given in Eq. (9) monotonically increases with At, 
therefore the derivative of the overlap function must be negative at the current 
time t (i.e., Vij -Tij < 0) for a interaction to take place now or at any time in the 
future. This condition is also the requirement that the particle must be before 
the root in the overlap function's derivative in its trajectory (marked by the tri- 
angle symbol in Fig. 3a). This condition also implicitly contains the condition 
that the particles have a non-zero relative velocity (|vy| > 0). Only the earliest 
root of the overlap function will generate interactions, as only the earliest root 
of the overlap function has a decreasing overlap function (marked with a filled 
circle in Fig. 3a). If the particles are already in contact or in an invalid state 
(rfj — afj < 0), then an interaction will take place immediately at time t as 
the first condition assures that f + (t) < (the particle pair is in the section of 
the trajectory marked with a solid line in Fig. 3a). If there is no overlap, the 
smallest real root of Eq. (8) is the time of the interaction. The calculation of 
the earliest root suffers from a computational problem known as catastrophic 
cancellation. In such a case, the alternate form of the quadratic equation must 
be used[19, 26]. If there are no real roots of the overlap function, the particles 
will miss each other, as illustrated in Fig. 3b. When the roots of Eq. (8) are a de- 
generate pair, no interaction occurs as the overlap function never changes sign. 
This completes the event detection for / + and a summary of the steps is listed 
in Algorithm 2. Unsurprisingly, this algorithm is similar to previously published 
algorithms [1, 11, 19] for detecting hard-sphere collisions; however, it differs by 
the addition of the second if-statement on line 3 of Algorithm 2. Crucially, it 
is this second if-statement which stabilises the event-detection algorithm with 
respect to invalid/overlapping states. 

1 Require = 0; 

2 if ■ r,ij > then return co; 

3 if tja — erf,; < then return t: 

4 if [(v^ ■ Tij) 2 - vf^rfj - afj)] < then return oo; 

s At <— (r?- - 4) / (-vn ■ nj + yWn,) 2 -^(4-4)) ; 
6 return t + At; 

Algorithm 2: The stable EDPD algorithm for event detection for the positive 
overlap function f£, =0 . This algorithm returns the time of the next transition 
across the positive discontinuity at cry for particle i and j. This algorithm is 
used to detect Core and Capture events in square-well models. 

The event-detection algorithm for interactions on the negative side of the dis- 
continuity, as illustrated in Fig. 3c-d, is now considered. This event-detection 
rule is used for Bounce/Release events in the square- well model and corresponds 
to the detection of the roots of the negative overlap function, /~ =0 , given in 
Eq. (8). The derivative of the overlap function is a monotonically decreasing 
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function with respect to At which implies there will always be an interaction 
provided vfj > 0, which forms the first interaction condition. For the root con- 
figuration given in Fig. 3c, it is again possible to solve for the time to interaction 
using the quadratic formula; however, if the argument of the square root is neg- 
ative, the particles do not cross the discontinuity at any point in their trajectory 
and we have the case illustrated in Fig. 3d. Even though the discontinuity is 
not crossed and the particle pair remains in an invalid state, interactions will 
still occur when the overlap function's derivative becomes negative. If the over- 
lap function's derivative is currently positive, an interaction must be scheduled 
at the time of the derivative's root (marked by a filled triangle in Fig. 3d). 
This corresponds to Line 10 of Algorithm 1 and is an illustration of the pre- 
viously discussed deficiency of the general algorithm. This case also occurs if 
the roots of the overlap function are degenerate. It is easy to see that in this 
case, any interaction scheduled at the time of the root in the overlap function's 
derivative would effectively be a "glancing" interaction and will not cause the 
overlap function to increase afterwards. Regardless of how this deficiency is 
handled, the stable algorithm requires that a interaction occurs immediately 
if f a . j=Q {t) < 0, or at a time t + At when the overlap function changes sign 

{f a . j=0 (t + At) = 0) to prevent the overlap function from decreasing while in 
an invalid state. In summary, the stable event-detection algorithm for f~ =0 
is listed in Algorithm 3. This algorithm is quite different to other published 
implementations [1] as it includes a case to handle the interactions of Fig. 3d on 
Line 8 and it also uses of a numerically stable form of the quadratic equation. 



1 Require a%j = 0: 

2 if vfj = then return oo; 

3 if ( Vij ■ rij f + v?. (afj - r?-) > then 



q < — •r i j + copysign | 

At <— max (0, (o£ - r?-) /q); 

return t + At; 

7 else 

8 I return At = max(0, — (v^j ■ fij)/vfj) ; 

9 end 

Algorithm 3: The stable EDPD algorithm for event detection for the negative 
overlap function /~ ._ . This algorithm returns the time of the next transition 
across the negative discontinuity at cry for particle i and j, It is used to detect 
Bounce/Release events in square well models. The function copysign(x, y) 
is a standard function which returns a value with the absolute value of the 
x expression and the sign of the y expression. The more common sign(y) 
function cannot be used as it returns zero for y = 0. 



Now that the optimal stable event-detection algorithms for particles with 
zero relative-acceleration have been defined, it is time to treat the more complex 
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Figure 4: The five classes of event-generating trajectories for non-zero relative accelerations 
(dij 7^ 0) of /+ . Only the configurations resulting in events have been drawn for brevity. The 
symbols and lines are described in the caption of Fig. 3. 



case of non-zero relative accelerations. This will require a strategy to solve for 
the roots of the quartic overlap functions. 

2.2. Non-zero Relative Acceleration (a^ ^0) 

It is desirable to have a stable event-detection algorithm for particle pairs 
with a non-zero relative acceleration, a,,- ^ 0, including both the positive f + 
and negative /" overlap functions. The positive overlap function is needed for 
particle systems which include gravity and use immobile spheres to form complex 
boundaries of the system[22] . The negative overlap function is required if square- 
bonds are to be used as constraints in rigid-body simulations or as bounding 
spheres for more complex shapes [24]. It is also technically interesting to derive 
these algorithms as they provide a non-trivial example for the application of the 
stable EDPD algorithm proposed in this work. 

As discussed before, the roots of quartic equations, such as the overlap func- 
tions in Eq. (6), can be solved using radicals; however, in general there is no 
strategy which is guaranteed to yield all roots[25]. Fortunately, such an algo- 
rithm is available for cubic equations [26] and so the roots of the time derivative 
of the overlap function (Eq. (7)) may be rapidly calculated. A stable event- 
detection routine can then be generated by bracketing roots in the overlap func- 
tion using the roots of it's derivative. It is now demonstrated that this strategy 
is sufficient to derive the interaction times of the positive overlap function. 

Considering only the trajectories which result in an interaction for the posi- 
tive overlap function, /+ , there are only five possible configurations of the roots 
and these are illustrated in Fig. 4. From the diagram it is obvious that there 
are only two possible regions in which interactions may occur. Either interac- 
tions occur between the current time and the first root in the overlap functions 
derivative, or between the second and third root (if they exist). The general 
strategy for finding interactions starts with the calculation of the roots of f + 
(Eq. (7)) using the algorithm for cubic equations [26]. There is either one root 
or three occurring at times t + At\, t + At2, and t + At^. Any roots which don't 
exist are assigned a value of positive infinity and the roots are then sorted so 
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1 Require Ojj ^ 0; 

2 Require n roots ; /* Count of the solutions to f~(t + At) = */ 

3 Require [A*i, Ai 2 , At 3 ] ; /* Sorted solutions to /~(i + Ai) = */ 

4 if Aii > and f+(t + Ah) < then 

5 if f+(t) < then 

6 return t; 

7 else 

8 return Biscct(/ + , f, t + Ati); 

9 end 

10 end 

11 if n roots = 3 and At 3 > and f + (t + At 3 ) < then 



At m i n = max(0, At 2 ); 
if f + (t + At mm ) < then 

return t + At m i n 
else 

return Bisect(/+, t + At mini t + Ai 3 ); 
end 



12 
13 
14 
15 
16 
17 

is end 

19 return oo; 

Algorithm 4: The stable EDPD algorithm for event detection for the positive 
overlap function / + . This algorithm returns the time of the next transition 
across the positive discontinuity at cr^ for particle i and j with non-zero rela- 
tive acceleration (a^ ^ 0). The function Bisect(/, ti, t 2 ) bisects for a root in 
the function /, where the sign change occurs in the interval \t\, ti)- 
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Figure 5: The seven classes of event-generating trajectories for non-zero relative accelerations 
(aij 0) of / _ . The symbols and lines are described in the caption of Fig. 3. 

that Aii < Ai2 < A*3. Provided the first root of the derivative is in the future 
(Aii > 0) and inside the discontinuity (/(< + Ati) < 0), the first region of inter- 
action is bracketed between the current time and t + At±. The algorithm must 
test if the particles have already passed the discontinuity (f(t) < 0) and issue 
an immediate interaction if so. These conditions are sufficient to detect the first 
region of interaction. The second region of interaction only exists if there are 
three roots of the overlap functions derivative, and the last root lies behind the 
discontinuity (f(t + At 3) < 0). Provided this root is in the future (At^ > 0), the 
interaction time is bracketed between i + A<2 and t + At^] however, if the second 
root of the overlap functions derivative is in the past (A<2 < 0), the current time 
should be used as a lower bound. If there is no sign change in this interval, the 
particle is already in the interaction region and an interaction should be issued 
immediately. This completes the description of the event-detection algorithm 
for the positive overlap function, / + , which is summarised in Algorithm 4. 

The final event-detection algorithm to consider is for the negative overlap 
function, f~. As the relative acceleration is non-zero, all trajectories will gener- 
ate interactions and they are illustrated in Fig. 5. As with the previous algorithm 
for the positive overlap function, the algorithm begins by calculating the time- 
ordered solutions to f~(t + At) = 0. There are again two distinct regions of the 
trajectory which may generate interactions. If / has three roots, one possible 
interaction region lies in the interval [A<i, A^) provided a sign change in the 
overlap function occurs (f(t + A^) < 0). If this region is in the past (A*2 < 0) 
or if there is only one root then the next interaction will occur in the interval 
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1 Require ^ 0; 

2 Require n roots ; /* Count of the solutions to f~(t + At) = */ 

3 Require [At x , Ai 2 , Ai 3 ] ; /* Sorted solutions to f~(t + At) = */ 
/* Catch the earliest interaction region of Figs. 5c — f */ 

4 if n roots = 3 and At2 > and f~(t + Ai 2 ) < then 

5 Ai min = max(0, Aii); 

6 if + Atmin) < then return t + At m m, 

7 return Bisect ( / ~~ , t + At m i n , t + At2 ) ; 

8 end 

9 Atust = Ai„ roots ; 

/* Catch the later interaction region of Figs. 5d,f ,g */ 

10 if f~(t + Atiast) < then return t + max(0, Ati as t); 

/* Handle all other cases */ 

n if Atiast < and f~(t) < then return t; 

12 return FindLastRoot(/~, i + max(0, At; as t)); 
Algorithm 5: The stable EDPD algorithm for event detection for the negative 
overlap function / _ . This algorithm returns the time of the next transition 
across the negative discontinuity at aij for particle i and j with non-zero rel- 
ative acceleration (a,j ^ 0). The Bisect function is described in Algorithm 4. 
The FindLastRoot(/, to) function searches for the sign change of the function 
f{t), from positive to negative value, occurring at some value after the initial 
guess time to- Only one root appears in this interval. 



[At„ roots , oo), where At nrootB corresponds to the last root of /. Unfortunately, 
the infinite upper limit of this search interval complicates the calculation of the 
interaction time as bisection is no longer possible. The procedure described so 
far and summarised in Algorithm 5 has reduced the detection for negative over- 
lap function events to the search for a FindLastRoot function. This function 
must find the single root of the overlap function which is guaranteed to occur 
in an interval with a finite lower bound. 

Although Householder root-finding methods such as the Newton-Raphson 
algorithm may be used in this case, such approaches can become unstable 
when the current time is near a root in the overlap function's derivatives (e.g., 
Ati — ► 0). There are more stable techniques available such as the homotopy 
methods [15] or those relying on the calculable bounds of the overlap functions 
derivatives in the search interval[14, 27]; however, as the overlap function has 
only a single root in the interval a much simpler approach is possible. 

Provided a reasonable estimate for the location of the root can be calculated, 
it can be rapidly bracketed using an exponentially growing search algorithm [26]. 
In the simplest case, where that the relative velocity and position are close to 
zero at the time of the final root in the overlap functions derivative (i+ Ai„ roots ), 
the negative overlap function given in Eq. (6) becomes 

At 4 

f-(t + At nroots + At) a 4 _ _ a 2._ 
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In this approximation, the latest root of / occurs at a time t + At nrootB + At e , 
where At e is given by the following expression 




This provides a rough order of magnitude estimate for the location of the final 
root of the overlap function. If the overlap function has changed sign at the time 
of this estimate (/~ (t + At nrootB + At e ) < 0), the root has been bracketed in the 
range [t+At nroots , t+At nroots +At e ). If not, the estimate of At e is doubled until 
a sign change is detected. An implementation of this FindLastRoot function 
using this approach is given in Algorithm 6. 



1 Require a. L j ^ 0; 

2 Require to ; /* The lower bound of the interval to search */ 

3 At e = (4 4A4) 1/4 ; 

4 while f-(t Q + At e ) > do 
s to = t + At e ; 

6 At e = 2At e ; 

7 end 

8 return Biscct(/~, to, t + At e ); 

Algorithm 6: An implementation of the FindLastRoot (/ _ , to) function, 
which searches for the single sign change of the negative overlap function 
f~(t), from positive to negative value, occurring at some time after the initial 
guess time to. The Bisect function is described in Algorithm 4. 



This completes the description of the stable event-detection routines for 
spherical potentials. Algorithm 2 has already been validated in high-precision 
simulations of the hard sphere gas [28]. A series of simulations of the square- 
well fluid have also been performed using Algorithm 3 [8]; therefore only the 
new Algorithms 4 and 5 are untested. Two numerically sensitive systems are 
simulated to test these new algorithms. The first simulation performed is of the 
gravity-driven flow of a square-well fluid from a hopper into a container (see 
Fig. 6a). The hopper, chute, and container are meshed using static dissipa- 
tive hard-sphere particles, which uses Algorithm 4 to detect collisions with the 
fluid particles. To test Algorithm 5, a particularly sensitive test case of a pen- 
dulum is simulated (see Fig. 6b). The pendulum chain consists of dissipativc 
hard-sphere particles which are bonded together using an infinite square-well 
interaction. The top of the chain is bonded using an infinite square-well to a 
static hard-sphere particle. Algorithm 5 is extensively tested in this case as it is 
used to implement the maximum separation constraint between the static parti- 
cle and the first chain particle. Both of these simulations are dissipative and the 
TC model[23] is used to avoid inelastic collapse. Each simulation was run over 
the full range of dissipative parameters, for millions of events, and no numer- 
ical difficulties were detected, validating the stability of these routines. These 
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a) b) 

Figure 6: Renderings of the two test systems for the stable EDPD algorithms: a) hopper flow 
of a square-well fluid and b) a pendulum. 

simulations help to highlight the complexity which may be captured using the 
analytical dynamics of EDPD and the event-detection algorithms proposed in 
this work. 

3. Conclusions 

A general and stable approach to event-detection in event-driven particle 
dynamics is proposed. It relies on treating interactions as stabilising events, 
and using these to treat invalid states which are not improving. Effectively, this 
defines the dynamics of invalid states to minimise their duration and their effect 
on the trajectory of the system. This leads to improved expressions for event 
detection, even for well-established discontinuous particle models. Algorithm 2 
and 3 are direct replacements for the common expressions used in hard-sphere, 
square-well, and stepped particle models. New expressions are also derived for 
complex particle systems involving varying particle accelerations, establishing 
the generality of the procedure and extending event-driven particle dynamics 
into these systems. Algorithm 4 may be used to implement hard-sphere meshes 
used for boundaries (see Fig. 6a) and Algorithm 5 allows distance constraints 
to be implemented (see Fig. 6b) in simulations. These results help to establish 
EDPD as an analytical technique for rigid-body simulations of spherical particle 
systems. 

Although this work focuses on spherical potentials, extensions to non-spherical 
systems may follow along the well-established lines of ellipsoidal particles [15] or 
terraced potentials [16]. The new expressions presented in this paper may be 
used to provide bounding-sphere neighbourhoods [24] which can be used as a 
search interval for numerical event-detection of complex particle shapes, such as 
spherocylinders or composite particles [19]. The stable event-detection algorithm 
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also readily generalises to more complex boundary objects, such as cylinders or 
planes; however, if the algorithm is extended to triangle meshes it will allow the 
simulation of complex geometries in biological processes [29] or the rapid design 
and import of boundaries from CAD programs for granular simulations [19]. 

All methods developed in this paper have been implemented in the open- 
source DynamO simulation package [30]. 
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